# ============================================================
# Thomas König & Stefan Eschenwecker
# The European Court of Justice and legal European integration
# 
# This script produces Figure 5 in the main text.
# ============================================================

# load packages
library(rstan)
library(brms)
library(tidyverse)
library(patchwork)

# plot options
theme_set(theme_bw(base_size = 15))

# load data
CourtData <- read.csv("Court/Data/PreparedCourtData.csv") %>%
  filter(all_not_applic == 0)

# load model
RegModel <- readRDS("Court/Models/UnequalPrecision.rds")


# QoI - Variance Ratios

# set up scenarios
Scenarios <- list(
  Observed = list(supra_signal_int = CourtData$supra_signal_int,
                  ms_signal_int = CourtData$ms_signal_int_std),
  EUWeakInt = list(supra_signal_int = "Weak Intensity"),
  MSMeanInt = list(ms_signal_int_std = 0),
  MSNoEUStrong = list(ms_signal_int_std = min(CourtData$ms_signal_int_std),
                      supra_signal_int = "Strong Intensity"),
  MSStrongEUStrong = list(ms_signal_int_std = 2,
                        supra_signal_int = "Strong Intensity"),
  MSStrongEUNo = list(ms_signal_int_std = 2,
                      supra_signal_int = "No Intensity")
)

ScenarioData <- map(names(Scenarios), ~ {
  name <- .x
  params <- Scenarios[[.x]]
  new_data <- CourtData
  for (col in names(params)) {
    new_data[[col]] <- params[[col]]
  }
  new_data$comparison_base <- name
  return(new_data)
}) %>%
  set_names(names(Scenarios))

# generate fitted values over scenarios
# (Hint: takes some time to run)
set.seed(1808)
PredictedDispersion <- map(ScenarioData, ~ {
  fitted(RegModel, re_formula = NULL, newdata = .x,
         summary = F, dpar = "disc", scale = "response")
})

# transform fitted values from log scale to variances
TransformedVariance <- map(PredictedDispersion, ~ (1 / exp(.x))^2)

# specify comparisons
comparisons <- list(
  c("Observed", "EUWeakInt"),
  c("Observed", "MSMeanInt"),
  c("MSStrongEUStrong", "MSStrongEUNo"),
  c("MSStrongEUStrong", "MSNoEUStrong")
)

# calculate variance ratios
VarRatios <- map(comparisons, ~ {
  var1 <- TransformedVariance[[.x[[1]]]]
  var2 <- TransformedVariance[[.x[[2]]]]
  return(var1 / var2)
}) %>%
  set_names(map_chr(comparisons, ~ paste(.x[[1]], "vs.", .x[[2]])))


# summarize by collapsing over draws and store in single df
FinalDF <- map_df(names(VarRatios), ~ {
  ratio_matrix <- VarRatios[[.x]]
  mean_val <- apply(ratio_matrix, 2, mean)
  low_ci <- apply(ratio_matrix, 2, quantile, probs = 0.025) 
  up_ci <- apply(ratio_matrix, 2, quantile, probs = 0.975) 
  tibble(
    mean = mean_val,
    low_ci = low_ci,
    up_ci = up_ci,
    comparison = factor(.x)
  )
})


# create plot of variance ratios over CDFs
VarPlot <- ggplot(FinalDF) +
  stat_ecdf(aes(x = mean), geom = "step", linewidth = .8) +
  stat_ecdf(aes(x = low_ci), geom = "step",
            linewidth = .8, linetype = "dotted") +
  stat_ecdf(aes(x = up_ci), geom = "step",
            linewidth = .8, linetype = "dotted") +
  geom_vline(xintercept = 1, linetype = "dashed",
             color = "grey60", linewidth = 1) +
  coord_flip() +
  facet_wrap(~comparison, ncol = 2, 
             labeller = labeller(comparison = c(
               "Observed vs. EUWeakInt" = "Weak EU Strength Baseline",
               "Observed vs. MSMeanInt" = "Mean MS Strength Baseline",
               "MSStrongEUStrong vs. MSStrongEUNo" = "Additional Effect of Strong EU Signal",
               "MSStrongEUStrong vs. MSNoEUStrong" = "Additional Effect of Strong MS Signal"
             ))) +
  xlab("Variance Ratio") +
  ylab("Posterior Cumulative Probability") +
  xlim(c(0, 8)) +
  scale_y_continuous(labels = scales::percent_format(), 
                     limits = c(0.001, 0.999))

# save plot
ggsave(plot = VarPlot, device = "pdf",
       filename = "VarianceRatios.pdf",
       path = "Court/Plots/", height = 9, width = 9)
